pro initialH
;;writed by wuxk 2019-1-18
print,'for loop must be used with'
print,'pro or function, carefully!'
;;needed to compile

g=file_import("data/circular.grd.hl2a22535.831ms.68x256.0.65to1.1.nc")
NX=g.NX
pp=fltarr(NX)
ni=fltarr(NX)
ti=fltarr(NX)
te=fltarr(NX)
pei=fltarr(NX)

Diff=fltarr(NX)
chii=fltarr(NX)
chie=fltarr(NX)
unit_psi=g.rmag*g.rmag*g.bmag
density_unit = 1.e19                          ; 1/m^3
lbar=g.rmag
ee = 1.6022e-19                               ; elementary charge 
Ni_x = 1
Te_x = 10                                     ; eV


jy=153                                        ;fueling path

rxy=g.Rxy[*,jy]/lbar
bpxy=g.Bpxy[*,jy]/g.bmag

psi=(g.psixy[*,jy]-g.psi_axis)/(g.psi_bndry-g.psi_axis)
xreal=g.rxy[*,jy]/lbar
print,'xreal(0),xreal(NX-1) = ',xreal(0),',',xreal(NX-1)
;xreal=g.psixy[*,jy]/unit_psi
;x=(g.rxy[*,jy]-min(g.rxy))/(max(g.rxy)-min(g.rxy))

p_core = 950.            ; pmin PW Xi in eV
p_edge = 2.              ; pmin - pbottom in PW Xi in eV
p_gradient = 0.01        ; parameter for gradient at pedestal top
p_enlarge =3.8           ; parameter for enlarge of core value
delta_ped=0.062          ; width of pedestal
x0_p = 0.90              ; center of pedestal relative to normalzied psi

;*** same  as dens5   H4
;ni_core = 0.075
;ni_edge = 0.048
;ti_core = 4.5
;ti_edge = 0.38
;te_core = 4.5
;te_edge = 0.38
  
;*** relative to dens6  H6
ni_core = 0.335          
ni_edge = 0.0111         
ti_core = 202.                 ;eV
ti_edge = 5.                    ;eV
te_core = ti_core               ;eV
te_edge = ti_edge               ;eV

Diff_xin = 1.          ; m^2/s
;Ga_xin=-Diff_xin*aveYg11J_core*dNidx_xin_au
;Qe_xin=-Diff_xin*aveYg11J_core*dTedx_xin_au
;Qi_xin=-Diff_xin*aveYg11J_core*dTidx_xin_au

print,'debug>>>>for loop is going to start!'
print,'debug>>>>psi(0), psi(NX-1)='+string(psi(0))+','+string(psi(NX-1))
;for jx=0,NX-1,1 DO print,string(psi(jx))
for jx=0,NX-1,1  do begin
     ;print,'jx =',jx
     xprm = (x0_p - psi(jx))/delta_ped
     Hmode_tanh =(exp(-xprm)+(1.+p_gradient*xprm)*exp(xprm)*p_enlarge)/(exp(-xprm)+exp(xprm))
     pp[jx]= p_edge +p_core* (Hmode_tanh-1.) 
     if (jx eq 0) then print,'       jx      psi[jx]      xprm       pp[jx]'
     print,jx,psi[jx],xprm,pp[jx]

     ni[jx]= ni_edge +ni_core* (Hmode_tanh-1.) 
     ti[jx]= ti_edge +ti_core* (Hmode_tanh-1.) 
     te[jx]= te_edge +te_core* (Hmode_tanh-1.) 
endfor

pei=ni*(ti+te)
pei = pei*Ni_x*density_unit*ee*Te_x             ; SI unit Pascals

print,'Variable', '          Min',    '           Max'
print,'ni  : ', min(ni),'---', max(ni)
print,'ti  : ', min(ti),'---', max(ti)
print,'te  : ', min(te),'---', max(te)
print,'pei : ', min(pei),'---', max(pei)


dppdx=deriv(xreal,pei)
;print,dppdx
dnidx=deriv(xreal,ni)
dtidx=deriv(xreal,ti)
dtedx=deriv(xreal,te)
;print,'dnidx-dtidx=',dnidx-dtidx

print,'Gradient', '    Core',    '      Edge'
print,'ni : ', dnidx[0],'---', dnidx[NX-1]
print,'ti : ', dtidx[0],'---', dtidx[NX-1]
print,'te : ', dtedx[0],'---', dtedx[NX-1]



flux_xin_ni = -Diff_xin*dnidx[0]*rxy[0]*bpxy[0]
flux_xin_ti = -Diff_xin*dtidx[0]*rxy[0]*bpxy[0]
flux_xin_te = -Diff_xin*dtedx[0]*rxy[0]*bpxy[0]

for jx = 1,NX-1  do begin

Diff[jx] = -flux_xin_ni/rxy[jx]/bpxy[jx]/dnidx[jx]
chii[jx] = -flux_xin_ti/rxy[jx]/bpxy[jx]/dtidx[jx]
chie[jx] = -flux_xin_te/rxy[jx]/bpxy[jx]/dtedx[jx]

;if Diff[jx] GT 10. then BEGIN
 ;Diff[jx]=10.
;ENDIF

ENDFOR

Diff[0]=Diff_xin
chii[0]=Diff_xin
chie[0]=Diff_xin

print,'Diffusion', '        Min',    '          Max'
print,'Diff  : ', min(Diff),'---', max(Diff)
print,'chi_i : ', min(chii),'---', max(chii)
print,'chi_e : ', min(chie),'---', max(chie)

ddxdiff=deriv(xreal,Diff)
ddxchii=deriv(xreal,chii)
ddxchie=deriv(xreal,chie)

print,'Gradient', '        Core',    '          Edge'
print,'Diff_ni : ', ddxdiff[0],'---', ddxdiff[NX-1]
print,'chii    : ', ddxchii[0],'---', ddxchii[NX-1]
print,'chie    : ', ddxchie[0],'---', ddxchie[NX-1]

safe_colors, /first


!p.multi=[0,2,3,0,0]
!x.omargin=[45,4]
window,1,xsize=900,ysize=700,xpos=0,ypos=74
plot,psi,pp,/xst,chars=2,charthick=1,thick=2,title='P fitting[Pa]'
plot,psi,pei,/xst, title='P[pa]',chars=2,charthick=1,thick=2

plot,psi,Diff,/xst,title='Diffusion Coef[m^2/s]',chars=2,charthick=1,thick=2
plot,psi,chii,/xst,title='chii Coef[m^2/s]',chars=2,charthick=1,thick=2

plot,psi,ni,/xst,chars=2,charthick=1,thick=2,title='Ni[N0]'
plot,psi,ti,/xst,chars=2,charthick=1,thick=2,title='Ti[eV]'
textx=0.01
texty=0.9
dy=0.05

xyouts, textx,texty,/normal,'diffusion_coef_Hmode0',charsize=2.5
xyouts, textx,texty-dy,/normal,' ='+string(Diff_xin),charsize=2.5

xyouts, textx,texty-2*dy,/normal,'Te_core='+string(Te_core),charsize=2
xyouts, textx,texty-3*dy,/normal,'Te_edge='+string(Te_edge),charsize=2
xyouts, textx,texty-4*dy,/normal,'Ni_core='+string(Ni_core),charsize=2
xyouts, textx,texty-5*dy,/normal,'Ni_edge='+string(Ni_edge),charsize=2

xyouts, textx,texty-7*dy,/normal,'p_position='+string(x0_p),charsize=2
xyouts, textx,texty-8*dy,/normal,'p_width   ='+string(delta_ped),charsize=2
xyouts, textx,texty-9*dy,/normal,'p_gradient='+string(p_gradient),charsize=2
xyouts, textx,texty-10*dy,/normal,'p_enlarge ='+string(p_enlarge),charsize=2

xyouts, textx,texty-12*dy,/normal,'flux_ni_xin ='+string(flux_xin_ni),charsize=2
xyouts, textx,texty-13*dy,/normal,'flux_ti_xin ='+string(flux_xin_ti),charsize=2
xyouts, textx,texty-14*dy,/normal,'flux_te_xin ='+string(flux_xin_te),charsize=2
xyouts, textx,texty-15*dy,/normal,'min diff ='+string(min(Diff)),charsize=2
xyouts, textx,texty-16*dy,/normal,'min chii ='+string(min(chii)),charsize=2

;window,2
;plot,psi,Diff,/xst,title='Diffusion Coef',chars=2,charthick=2,thick=2

;window,3
;plot,x,chii,/xst,title='chii Coef',chars=2,charthick=2,thick=2

;window,4
;plot,x,chie,/xst,title='chie Coef',chars=2,charthick=2,thick=2

;window,5
;plot,x,ni,/xst,chars=2,charthick=2,thick=2,title='Ni'

;window,6
;plot,x,ti,/xst,chars=2,charthick=2,thick=2,title='Ti'

;window,7
;plot,x,te,/xst,chars=2,charthick=2,thick=2,title='Te'

;window,8
;plot,x,pei,/xst, title='P[pa]',chars=2,charthick=2,thick=2

end


